% Comment out the next line if the program is used in 'simulation mode' without SMM annealing loop.
function [criterion simmom IK_use HL_use K_use L_use Y_use V_use A_use Optk Optl lKmin lLmin lKmax lLmax dK_use dL_use] = simulation(params, nK, nL, nA, nF, extraA, Lmode, FC, mode, momlist)

warning off all

% ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------    
%                  List of arguments, simulation.m
% ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------  
%  params = Parameters of the model (adj. costs and AR(1) process).
%  nK     = Number of elements in Capital vectors.
%  nL     = Number of elements in Labor vectors.
%  nA     = Number of elements in finite state support for the AR(1).
%  nF     = Number of possible plant types.
%  extraA = Number of extra elements in finite state support for interpolation after value function iteration.
%  Lmode  = Determines whether labor fixed costs are associated with gross or net flows.   
%  FC     = Determines whether fixed costs are considered scaled or purely fixed.
%  mode   = Parameter search ('search') or single simulation ('Sim').
%  momlist = Determines which moment vector is the basis for evaluation criterion.   
% ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------     

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %             Define parameters             % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

beta       = 0.95;                                                          % Yearly discount rate.
dK         = 0.1;                                                           % Capital depreciation rate (mean).
dL         = 0.1;                                                           % Labor separation rate (mean). 
pI         = 1.0;                                                           % Unit purchase price of capital.  
w1         = 1.00;                                                          % Parameter of the wage function: cost per worker per year.
etaK       = .25;                                                           % Capital output elasticity.
etaL       = .50;                                                           % Labor output elasticity.
mu         = 0.0;                                                           % Mean shock level in firm level shocks.
A_bar      = 1.0;                                                           % Median productivity level. 0.0755
bK         = params(1);                                                     % Convex cost parameter, capital.
bL         = params(2);                                                     % Convex cost parameter, labor.
alphaK     = params(3);                                                     % Fixed cost parameter, capital.
alphaL     = params(4);                                                     % Fixed cost parameter, labor.
alphaKL    = params(5);                                                     % Interrelation cost parameter.
pH         = 0;                                                             % Linear cost compone 0.2739nt to changing number of workers.
resaleloss = params(6);                                                     % Share of capital value lost when making negative investments.
rho        = params(7);                                                     % Shock persistence in AR(1) process for firm level shocks.
sigma      = params(8);                                                     % Standard dev. in AR(1) process for firm level shocks.
sigdK      = params(9);                                                     % Spread in distribution of firm types (cap. dep.).
sigdL      = params(10);                                                    % Spread in distribution of firm types (lab. sep.).
ro         = params(11);                                                    % Persistency in firm level heterogeneity (should be very close or equal to one).  
me(1)      = params(12);                                                    % Variance of capital measurement error.
me(2)      = params(13);                                                    % Variance of labor measurement error.

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %       Check valid function arguments      % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

if strcmp(Lmode,'LaborNet') + strcmp(Lmode,'LaborGross') == 0;
 fprintf('Invalid mode for labor fixed costs. Specify LaborNet or LaborGross \n')
elseif strcmp(FC,'Fixed') + strcmp(FC,'Scaled') == 0;
 fprintf('Invalid mode for fixed costs. Specify Fixed or Scaled \n')
elseif strcmp(mode,'Sim') + strcmp(mode,'Search') == 0;
 fprintf('Invalid execution mode. Specify Sim or Search \n')
else
  

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % %  Parameters for the state space and grid  % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

kinc    = 0.1/3;                                                            % Difference (in logs) between gridpoints in capital vector.
linc    = 0.1/3;                                                            % Difference (in logs) between gridpoints in labor vector.
k0      = -1.90;                                                            % Minimum (mean) level of log(capital).
l0      = -2.70;                                                            % Minimum (mean) level of log(labor).
aK      = round((0.0*sigma + 0.0*rho^2)/exp(0.5*bK + 0.2*bL)^0.5);          % Difference btw mean and minimum shock type realisation wrt capital vector values.
aL      = round((0.0*sigma + 0.0*rho^2)/exp(0.2*bK + 0.5*bL)^0.5);          % Difference btw mean and minimum depreciation firm type realisation wrt capital vector values.
fK      = floor((0.0*sigdK + 0.0*sigdL)*1e4*ro);                            % Difference btw mean and minimum shock type realisation wrt labor vector values.
fL      = floor((0.0*sigdK + 0.0*sigdL)*1e4*ro);                            % Difference btw mean and minimum dpreciation firm type realisation wrt labor vector values.
K0(1,1) = k0 - floor(nA/2)*kinc*aK + floor(nF/2)*kinc*fK;                   % Minumum level of log(capital) at lowest type shock and low type depreciation.
L0(1,1) = l0 - floor(nA/2)*linc*aL + floor(nF/2)*linc*fL;                   % Minimum level of log(labor) at lowest type shock and low type depreciation.
K0(2,1) = k0 + floor(nA/2)*kinc*aK + floor(nF/2)*kinc*fK;                   % Minimum level of log(capital) at highest type shock and lowest type depreciation.
L0(2,1) = l0 + floor(nA/2)*linc*aL + floor(nF/2)*linc*fL;                   % Minimum level of log(labor) at highest type shock and lowest type depreciation.
K0(1,2) = k0 - floor(nA/2)*kinc*aK - floor(nF/2)*kinc*fK;                   % Minimum level of log(capital) at lowest type shock and highest type depreciation.
L0(1,2) = l0 - floor(nA/2)*linc*aL - floor(nF/2)*linc*fL;                   % Minimum level of log(labor) at lowest type shock and highest type depreciation.
K0(2,2) = k0 + floor(nA/2)*kinc*aK - floor(nF/2)*kinc*fK;                   % Minimum level of log(capital) at highest type shock and highest type depreciation.
L0(2,2) = l0 + floor(nA/2)*linc*aL - floor(nF/2)*linc*fL;                   % Minimum level of log(labor) at highest type shock and highest type depreciation.

aKmin(1,:) = linspace(K0(1,1),K0(2,1),nA);                                  % Minimum levels of log(capital) as a vector corresponding to the finite state space support for the AR(1) process for the lowest type of depreciation.
aLmin(1,:) = linspace(L0(1,1),L0(2,1),nA);                                  % Minimum levels of log(labor) as a vector corresponding to the finite state space support for the AR(1) process for the lowest type of depreciation. 
aKmin(2,:) = linspace(K0(1,2),K0(2,2),nA);                                  % Minimum levels of log(capital) as a vector corresponding to the finite state space support for the AR(1) process for the highest type of depreciation
aLmin(2,:) = linspace(L0(1,2),L0(2,2),nA);                                  % Minimum levels of log(labor) as a vector corresponding to the finite state space support for the AR(1) process for the highest type of depreciation.
lKmin      = ndlinspace(aKmin(1,:),aKmin(2,:),nF);                          % Minimum levels of log(capital) as a matrix corresponding to productivity and firm types.
lLmin      = ndlinspace(aLmin(1,:),aLmin(2,:),nF);                          % Minimum levels of log(labor) as a matrix corresponding to productivity and firm types.
lKmax      = lKmin + kinc*(nK-1);                                           % Maximum level of log(capital).
lLmax      = lLmin + linc*(nL-1);                                           % Maximum level of log(labor).
extraK     = 2*(1+(nA-1)*aK+(nF-1)*fK);                                     % Number of extra nodes at which to extrapolate values outside of capital grid.
extraL     = 2*(1+(nA-1)*aL+(nF-1)*fL);                                     % Number of extra nodes at which to extrapolate values outside of labor grid. 
detail     = 4;                                                             % Number of extra nodes between each element on the grid in interpolated policy function.
K          = exp(ndlinspace(lKmin', lKmax', nK));                           % Capital matrix.
L          = exp(ndlinspace(lLmin', lLmax', nL));                           % Labor matrix.
K          = permute(reshape(K',nK,nA,nF),[2 1 3]);                         % Reshape capital matrix.
L          = permute(reshape(L',nK,nA,nF),[2 1 3]);                         % Reshape labor matrix.

% To allow adjustments to levels outside K, we create the matrices K_long and L_long (size: nA x (nK+extraK) and nA x (nL+extraL)).   
K_long = [permute(reshape(exp(ndlinspace(lKmin'-(extraK/2)*kinc,lKmin',extraK/2+1)),nA,nF,extraK/2+1),[1 3 2]),K(:,2:end-1,:),permute(reshape(exp(ndlinspace(lKmax',lKmax'+(extraK/2)*kinc,extraK/2+1)),nA,nF,extraK/2+1),[1 3 2])];
L_long = [permute(reshape(exp(ndlinspace(lLmin'-(extraL/2)*linc,lLmin',extraL/2+1)),nA,nF,extraL/2+1),[1 3 2]),L(:,2:end-1,:),permute(reshape(exp(ndlinspace(lLmax',lLmax'+(extraL/2)*linc,extraL/2+1)),nA,nF,extraL/2+1),[1 3 2])];

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % % % %  Parameters for the simulation  % % % % % % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

N      = 1794;                                                              % Number of firms.
T      = 110;                                                               % Number of periods.
Tlast  = 10;                                                                % Number of periods to keep.
Kappa  = 10;                                                                % Number of simulation rounds.
seed   = 231;                                                               % Set seed value to be equal for all loops in the annealing algorithm.

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % Load vector of empirical moments and covariance matrix from file. % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

load('empmom.mat', 'EY', 'VY');
empmom = EY;               
covar  = VY;                

% Call the routine that gives us AR process with a finite state Markov chain.
[Transfirm, ~, af]              = tauchen(nA,mu,rho,sigma);        
[~, Epsfirmlong, aflong] = tauchen(nA+extraA,mu,rho,sigma);

% If allowing for firm specific depreciation rates, create new vector of depr rates + transition matrix.
if nF > 1;
 if ro < 1.0;
  [~, Epsdepr, dK]   = tauchen(nF, dK, ro, sigdK);
  [Transtype, ~, dL] = tauchen(nF, dL, ro, sigdL);
 else
  [~, Epsdepr, dK]   = tauchen(nF, dK, ro-0.0001, sigdK);
  [~, ~, dL]         = tauchen(nF, dL, ro-0.0001, sigdL);
  Transtype = eye(nF);
 end
else
 Transtype = 1;
 Epsdepr   = 0;
end

% Generate lognormal productivity vector based on the finite states in af and am, and scale productivity by A_bar.
% -> Size of productivity matrix A is (# of states in firm level shocks) x (# of states in macro level shocks) x (# of states in time-invariant effects).   
A     = exp(af).*A_bar;
Along = exp(aflong).*A_bar;

% Matrices containing binary indicator for, nonzero investment rates and hirings.
i           = single(repmat(permute(K_long,[2 4 1 3]),[1 nK 1 1])./repmat(permute(K,[4 2 1 3]),[(nK+extraK) 1 1])-(1-repmat(permute(dK,[2 3 4 1]),[(nK+extraK) nK nA 1])));
h           = single(repmat(permute(L_long,[2 4 1 3]),[1 nL 1 1])./repmat(permute(L,[4 2 1 3]),[(nL+extraL) 1 1])-(1-repmat(permute(dL,[2 3 4 1]),[(nL+extraL) nL nA 1])));
i_nonzero   = abs(i) > 0.015;
i_neg       = i < 0.0;
if strcmp(Lmode,'LaborNet') == 1;
 dL_nonzero = (h > repmat(permute(dL,[2 3 4 1]),[(nL+extraL) nL nA 1]) + .015) + (h < repmat(permute(dL,[2 3 4 1]),[(nL+extraL) nL nA 1]) - .015);
else
 dL_nonzero = abs(h)>0.015;
end
% Create array of returns under dimensions {capital next, labor next, capital now, labor now, productivity, firm type}. 

if strcmp(FC,'Fixed') == 1;
 returns = single(repmat(permute(A,[2 3 4 5 1]),[(nK+extraK) (nL+extraL) nK nL 1 nF]).*repmat(permute(K,[4 5 2 6 1 3]),[(nK+extraK) (nL+extraL) 1 nL 1 1]).^etaK.*repmat(permute(L,[4 5 6 2 1 3]),[(nK+extraK) (nL+extraL) nK 1 1 1]).^etaL ...    % Production/sales value.
         - repmat(permute(repmat(permute(K,[4 2 1 3]),[(nK+extraK) 1 1 1]).*((1/2).*bK.*i.^2 + pI*i - resaleloss.*i_neg.*i),[1 5 2 6 3 4]),[1 (nL+extraL) 1 nL 1 1]) ...                                                                            % Convex and direct costs/revenues from aqcuiring/selling capital.
         - repmat(permute(alphaK.*i_nonzero,[1 5 2 6 3 4]),[1 (nL+extraL) 1 nL 1 1]) ...                                                                                                                                                             % Fixed costs of capital adjustments.
         - repmat(permute(repmat(permute(L,[4 2 1 3]),[(nL+extraL) 1 1 1]).*((1/2).*bL.*h.^2 + pH.*h),[5 1 6 2 3 4]),[(nK+extraK) 1 nK 1 1 1]) ...                                                                                                  % Convex, and linear costs of gross net labor adjustments.
         - repmat(permute(alphaL.*dL_nonzero,[5 1 6 2 3 4]),[(nK+extraK) 1 nK 1 1 1]) ...                                                                                                                                                            % Fixed costs of labor adjustments.
         - repmat(permute(w1.*L,[4 5 6 2 1 3]),[(nK+extraK) (nL+extraL) nK 1 1 1]) ...                                                                                                                                                              % Wage costs.
         - repmat(permute(alphaKL.*repmat(permute(i_nonzero, [1 5 2 6 3 4]),[1 (nL+extraL) 1 nL 1 1]).*repmat(permute(dL_nonzero,[5 1 6 2 3 4]),[(nK+extraK) 1 nK 1 1]),[1 2 3 4 5 6]),[1 1 1 1 1 1]));  % Interrelated costs/discount of simultaneous adjustments
else
 returns = single(repmat(permute(A,[2 3 4 5 1]),[(nK+extraK) (nL+extraL) nK nL 1 nF]).*repmat(permute(K,[4 5 2 6 1 3]),[(nK+extraK) (nL+extraL) 1 nL 1 1]).^etaK.*repmat(permute(L,[4 5 6 2 1 3]),[(nK+extraK) (nL+extraL) nK 1 1 1]).^etaL ...    % Production/sales value.
         - repmat(permute(repmat(permute(K,[4 2 1 3]),[(nK+extraK) 1 1 1]).*((1/2).*bK.*i.^2 + alphaK.*i_nonzero + pI*i - resaleloss.*i_neg.*i),[1 5 2 6 3 4]),[1 (nL+extraL) 1 nL 1 1]) ...                                                       % Convex and fixed costs and direct costs/revenues from aqcuiring/selling capital.
         - repmat(permute(repmat(permute(L,[4 2 1 3]),[(nL+extraL) 1 1 1]).*((1/2).*bL.*h.^2 + alphaL.*dL_nonzero + pH.*h),[5 1 6 2 3 4]),[(nK+extraK) 1 nK 1 1 1]) ...                                                                            % Convex, linear and fixed costs of gross labor adjustments.
         - repmat(permute(w1.*L,[4 5 6 2 1 3]),[(nK+extraK) (nL+extraL) nK 1 1 1]) ...                                                                                                                                                             % Wage costs.
         - repmat(permute(alphaKL.*repmat(permute(i_nonzero, [1 5 2 6 3 4]),[1 (nL+extraL) 1 nL 1 1]).*repmat(permute(dL_nonzero,[5 1 6 2 3 4]),[(nK+extraK) 1 nK 1 1]).*repmat(permute(K,[4 5 2 6 1 3]),[(nK+extraK) (nL+extraL) 1 nL 1 1]).^(1/2).*repmat(permute(L,[4 5 6 2 1 3]),[(nK+extraK) (nL+extraL) nK 1 1 1]).^(1/2),[1 2 3 4 5 6]),[1 1 1 1 1 1]));  % Interrelated costs/discount of simultaneous adjustments
end

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % %    Start the value function iteration procedure.    % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %    
    
V         = ones(nK,nL,nA,nF);                                              % Value function initial guess.
V_extrap  = zeros(nK+extraK,nL+extraL,nA,nF);                               % Zero matrix.    
VV_extrap = zeros(nK+2,nL+2,nA,nA,nF,nF);                                   % Zero matrix.
Maxiter = 150;                                                              % Maximum number of value function iterations before returning NaN if convergence criterion not obtained.
ii = 0;                                                                      % Number of iterations so far is 0.

while ii < Maxiter;                
  % Start each iteration by extrapolating on each end of the current guess for the value function. Simply use bilinear extrapolation. 
  for j=1:nF; 
   for f=1:nA;
    V_extrap(1:(extraK/2),1:(nL+extraL),f,j)                = [((K(f,2,j)-K_long(f,1:extraK/2,j))'*(L(f,2,j)-L_long(f,1:extraL/2,j)).*V(1,1,f,j) + (K_long(f,1:extraK/2,j)-K(f,1,j))'*(L(f,2,j)-L_long(f,1:extraL/2,j)).*V(2,1,f,j) + (K(f,2,j)-K_long(f,1:extraK/2,j))'*(L_long(f,1:extraL/2,j)-L(f,1,j)).*V(1,2,f,j) + (K_long(f,1:extraK/2,j)-K(f,1,j))'*(L_long(f,1:extraL/2,j)-L(f,1,j)).*V(2,2,f,j))./((K(f,2,j)-K(f,1,j)).*(L(f,2,j)-L(f,1,j)))    ,((K_long(f,1:extraK/2,j)-K(f,1,j))/(K(f,2,j)-K(f,1,j)))'*(V(2,:,f,j)-V(1,:,f,j)) + repmat(V(1,:,f,j),[extraK/2 1]),      ((K(f,2,j)-K_long(f,1:extraK/2,j))'*(L(f,end,j)-L_long(f,nL+extraL/2+1:nL+extraL,j)).*V(1,end-1,f,j) + (K_long(f,1:extraK/2,j)-K(f,1,j))'*(L(f,end,j)-L_long(f,nL+extraL/2+1:nL+extraL,j)).*V(2,end-1,f,j) + (K(f,2,j)-K_long(f,1:extraK/2,j))'*(L_long(f,nL+extraL/2+1:nL+extraL,j)-L(f,end-1,j)).*V(1,end,f,j) + (K_long(f,1:extraK/2,j)-K(f,1,j))'*(L_long(f,nL+extraL/2+1:nL+extraL,j)-L(f,end-1,j)).*V(2,end,f,j))./((K(f,2,j)-K(f,1,j)).*(L(f,end,j)-L(f,end-1,j)))];
    V_extrap((extraK/2+1):(nK+extraK/2),1:(nL+extraL),f,j)  = [(V(:,2,f,j)-V(:,1,f,j))*((L_long(f,1:extraL/2,j)-L(f,1,j))/(L(f,2,j)-L(f,1,j))) + repmat(V(:,1,f,j),[1 extraL/2]), V(:,:,f,j), (V(:,end,f,j)-V(:,end-1,f,j))*((L_long(f,nL+extraL/2+1:end,j)-L(f,end,j))/(L(f,end,j)-L(f,end-1,j))) + repmat(V(:,end,f,j),[1 extraL/2])];
    V_extrap((nK+extraK/2+1):(nK+extraK),1:(nL+extraL),f,j) = [((K(f,end,j)-K_long(f,nK+extraK/2+1:nK+extraK,j))'*(L(f,2,j)-L_long(f,1:extraL/2,j)).*V(end-1,1,f,j) + (K_long(f,nK+extraK/2+1:nK+extraK,j)-K(f,end-1,j))'*(L(f,2,j)-L_long(f,1:extraL/2,j)).*V(end,1,f,j) + (K(f,end,j)-K_long(f,nK+extraK/2+1:nK+extraK,j))'*(L_long(f,1:extraL/2,j)-L(f,1,j)).*V(end-1,2,f,j) + (K_long(f,nK+extraK/2+1:nK+extraK,j)-K(f,end-1,j))'*(L_long(f,1:extraL/2,j)-L(f,1,j)).*V(end,2,f,j))./((K(f,end,j)-K(f,end-1,j)).*(L(f,2,j)-L(f,1,j)))    ,((K_long(f,nK+extraK/2+1:end,j)-K(f,end,j))./(K(f,end,j)-K(f,end-1,j)))'*(V(end,:,f,j)-V(end-1,:,f,j)) + repmat(V(end,:,f,j),[extraK/2 1]),     ((K(f,end,j)-K_long(f,nK+extraK/2+1:nK+extraK,j))'*(L(f,end,j)-L_long(f,nL+extraL/2+1:nL+extraL,j)).*V(end-1,end-1,f,j) + (K_long(f,nK+extraK/2+1:nK+extraK,j)-K(f,end-1,j))'*(L(f,end,j)-L_long(f,nL+extraL/2+1:nL+extraL,j)).*V(end,end-1,f,j) + (K(f,end,j)-K_long(f,nK+extraK/2+1:nK+extraK,j))'*(L_long(f,nL+extraL/2+1:nL+extraL,j)-L(f,end-1,j)).*V(end-1,end,f,j) + (K_long(f,nK+extraK/2+1:nK+extraK,j)-K(f,end-1,j))'*(L_long(f,nL+extraL/2+1:nL+extraL,j)-L(f,end-1,j)).*V(end,end,f,j))./((K(f,end,j)-K(f,end-1,j)).*(L(f,end,j)-L(f,end-1,j)))];
   end
  end
 
  % Dimensions of expected value matrix are {capital next, labor next, expected firm level shock, expected firm type}.
  % We find the expected value by taking expectations over idiosyncratic firm level shocks and realization of firm type transition. Matlab function dot calculates vectorproduct in specified dimension.  
  for f = 1 : nF
   for ff = 1 : nF
    for a = 1 : nA
     for aa = 1 : nA
      VV_extrap(:,:,a,aa,f,ff) = permute(V_extrap(extraK/2+(a-aa)*aK+(f-ff)*fK:end-extraK/2+1+(a-aa)*aK+(f-ff)*fK,extraL/2+(a-aa)*aL+(f-ff)*fL:end-extraL/2+1+(a-aa)*aL+(f-ff)*fL,aa,ff),[1 2 5 3 6 4]);
     end
    end
   end
  end
  EV_extrap = squeeze(single(dot(repmat(permute(Transfirm,[3 4 1 2]),[(nK+2) (nL+2) 1 1 nF nF]),VV_extrap,4)));
  EV_extrap = squeeze(single(dot(repmat(permute(Transtype,[3 4 5 1 2]),[(nK+2) (nL+2) nA 1 1]),EV_extrap,5)));

  % Take new guess at the expected value for all possible choices = returns + value next period. Next command performs 3 operations in 1 go: 1)Use replicate and permute commands to expand EV_extrap matrix to same dimensions as returns matrix. 2)Set the temporary value function matrix to the value corresponding to optimal investment given today's level of capital. Start by maximizing over capital dimension, then labor. 3)Get rid of singleton dimensions in TV using squeeze function. After this, TV will have dimensionality {nK, nL, length(Prob)}.    
  TV = squeeze(squeeze(max(max(single(returns(extraK/2:nK+extraK/2+1,extraL/2:nL+extraL/2+1,:,:,:,:) + beta.*repmat(permute(EV_extrap,[1 2 6 7 3 4 5]),[1 1 nK nL 1 1 1])),[],1),[],2))); 
  % Evalueate the difference between current and previous value function guess (just use one of the productivity levels to evaluate).    
  % Update the value function guesstimates. 
  V=TV;
  % Update the iteration count + 1.
  ii = ii+1;
end

% After obtaining evaluation for chosen grid, get corresponding policy function.  
[TV Optk] = max(single(returns(extraK/2:nK+extraK/2+1,extraL/2:nL+extraL/2+1,:,:,:,:) + beta.*repmat(permute(EV_extrap,[1 2 6 7 3 4 5]),[1 1 nK nL 1 1])),[],1);
[~, Optl] = max(TV,[],2);
Tempoptk  = zeros(nK,nL,nA,nF);
for j=1:nF
 for f=1:nA
  for l=1:nL
   for k=1:nK
    Tempoptk(k,l,f,j) = Optk(1,Optl(1,1,k,l,f,j),k,l,f,j);
   end
  end
 end
end
Optk = Tempoptk;
Optl = squeeze(squeeze(Optl));
% Since Optk and Optl now correspond to gridpoints in K_long and L_long, respectively, we need to 'translate' these points into nodes on our initial K vector.
Optk = Optk - 1;
Optl = Optl - 1;

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % %               Start the simulations.                % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

% Continue the process of simulation only if the value function ended before j = MAXITER, and if optimal policies stay within grid.
% Otherwise, return criterion = NaN.
if (min(min(min(min(min(Optk)))))>0)*(max(max(max(max(max(Optk)))))<=nK)*(min(min(min(min(min(Optl)))))>0)*(max(max(max(max(max(Optl)))))<=nL) == 0;
 fprintf('----------------------------------------------------------- \n')
 fprintf(' Warning: Optimal input levels outside state space limits. \n')
 fprintf('----------------------------------------------------------- \n')
 fprintf('min(OptK) > 0: %1.0f  max(OptK) <= nK: %1.0f  min(OptL) > 0: %1.0f  max(OptL) <= nL: %1.0f \n',(min(min(min(min(min(Optk)))))>0),(max(max(max(max(max(Optk)))))<=nK),(min(min(min(min(min(Optl)))))>0),(max(max(max(max(max(Optl)))))<=nL))
 fprintf('----------------------------------------------------------- \n')
 fprintf('OptK(1,1,1,1)         = %3.0f, OptL(1,1,1,1)         = %3.0f \n', Optk(1,1,1,1), Optl(1,1,1,1))
 fprintf('OptK(1,1,1,end)       = %3.0f, OptL(1,1,1,end)       = %3.0f \n', Optk(1,1,1,end), Optl(1,1,1,end))
 fprintf('OptK(end,end,end,1)   = %3.0f, OptL(end,end,end,1)   = %3.0f \n', Optk(end,end,end,1), Optl(end,end,end,1))
 fprintf('OptK(end,end,end,end) = %3.0f, OptL(end,end,end,end) = %3.0f \n', Optk(end,end,end,end), Optl(end,end,end,end))
 fprintf('----------------------------------------------------------- \n')
end
if (min(Optk(1,1,1,:))>0)*(min(Optl(1,1,1,:))>0)*(max(Optk(end,end,end,:))<nK)*(max(Optl(end,end,end,:))<nL) == 0
 criterion = NaN('single'); simmom = NaN('single'); IK_use = NaN('single'); HL_use = NaN('single'); K_use = NaN('single'); L_use = NaN('single'); Y_use = NaN('single'); A_use = NaN('single'); V_use = NaN('single');
 fprintf('----------------------------------------------------------- \n')
 fprintf('Cannot compute simulations based on found policy functions. \n')
 fprintf('----------------------------------------------------------- \n')
else
    
% First, translate policy functions into investment and hiring rates.    
Kresh = repmat(permute(K,[2 4 1 3]),[1 nL 1 1]);
Lresh = repmat(permute(L,[4 2 1 3]),[nK 1 1 1]);
Kvect = reshape(Kresh,nK*nL*nA*nF,1,1,1);
Lvect = reshape(permute(Lresh,[2 1 3 4]),nK*nL*nA*nF,1,1,1);
Optik = Kvect(Optk+1)./Kresh - (1-repmat(permute(dK,[2 3 4 1]),[nK nL nA 1]))-kinc;
Opthl = Lvect(Optl+1)./Lresh - (1-repmat(permute(dL,[2 3 4 1]),[nK nL nA 1]))-linc;

% Create detailed policy functions by linear interpolation:
Kdet = exp(ndlinspace(lKmin', lKmax', nK+(nK-1)*detail));   
Ldet = exp(ndlinspace(lLmin', lLmax', nL+(nL-1)*detail));
Kdet = permute(reshape(Kdet',nK+(nK-1)*detail,nA,nF),[2 1 3]);
Ldet = permute(reshape(Ldet',nL+(nL-1)*detail,nA,nF),[2 1 3]);

Optkdetk     = zeros(nK+(nK-1)*detail,nL,nA,nF); Optldetk  = zeros(nK+(nK-1)*detail,nL,nA,nF);
Vdetk        = zeros(nK+(nK-1)*detail,nL,nA,nF); Optkdetkl = zeros(nK+(nK-1)*detail,nL+(nL-1)*detail,nA,nF); 
Optldetkl    = zeros(nK+(nK-1)*detail,nL+(nL-1)*detail,nA,nF); Vdetkl = zeros(nK+(nK-1)*detail,nL+(nL-1)*detail,nA,nF);
Optikdetk    = zeros(nK+(nK-1)*detail,nL,nA,nF); Optikdetkl = zeros(nK+(nK-1)*detail,nL+(nL-1)*detail,nA,nF);
Opthldetk    = zeros(nK+(nK-1)*detail,nL,nA,nF); Opthldetkl = zeros(nK+(nK-1)*detail,nL+(nL-1)*detail,nA,nF);

for j = 1:nF
 for i = 1:nA
  for l = 1:nL
   Optkdetk(:,l,i,j) = interp1(permute(K(i,:,j),[2 3 1]),Optk(:,l,i,j),permute(Kdet(i,:,j),[2 1 3]));
   Optldetk(:,l,i,j) = interp1(permute(K(i,:,j),[2 3 1]),Optl(:,l,i,j),permute(Kdet(i,:,j),[2 1 3]));
   Vdetk(:,l,i,j) = interp1(permute(K(i,:,j),[2 3 1]),V(:,l,i,j),permute(Kdet(i,:,j),[2 1 3]));
   Optikdetk(:,l,i,j) = interp1(permute(K(i,:,j),[2 3 1]),Optik(:,l,i,j),permute(Kdet(i,:,j),[2 1 3]));
   Opthldetk(:,l,i,j) = interp1(permute(K(i,:,j),[2 3 1]),Opthl(:,l,i,j),permute(Kdet(i,:,j),[2 1 3]));
  end
  for k = 1:length(Kdet)
   Optkdetkl(k,:,i,j) = interp1(permute(L(i,:,j),[2 3 1]),Optkdetk(k,:,i,j),permute(Ldet(i,:,j),[2 1 3]));
   Optldetkl(k,:,i,j) = interp1(permute(L(i,:,j),[2 3 1]),Optldetk(k,:,i,j),permute(Ldet(i,:,j),[2 1 3]));
   Vdetkl(k,:,i,j) = interp1(permute(L(i,:,j),[2 3 1]),Vdetk(k,:,i,j),permute(Ldet(i,:,j),[2 1 3]));
   Optikdetkl(k,:,i,j) = interp1(permute(L(i,:,j),[2 3 1]),Optikdetk(k,:,i,j),permute(Ldet(i,:,j),[2 1 3]));
   Opthldetkl(k,:,i,j) = interp1(permute(L(i,:,j),[2 3 1]),Opthldetk(k,:,i,j),permute(Ldet(i,:,j),[2 1 3]));
  end
 end
end


Vdetkla        = permute(interp1(af,permute(Vdetkl,[3 1 2 4]),aflong,'linear','extrap'),[2 3 1 4]);
IKdetkla       = permute(interp1(af,permute(Optikdetkl,[3 1 2 4]),aflong,'linear','extrap'),[2 3 1 4]);
HLdetkla       = permute(interp1(af,permute(Opthldetkl,[3 1 2 4]),aflong,'linear','extrap'),[2 3 1 4]);
Kdeta          = interp1(af,Kdet,aflong,'linear','extrap');
Ldeta          = interp1(af,Ldet,aflong,'linear','extrap');
for i = 1 : nA + extraA 
 IKdetkla(:,:,i,:)   = IKdetkla(:,:,i,:)*(Along(i)>A(1))*(Along(i)<A(end)) + Optikdetkl(:,:,1,:)*(Along(i)<A(1)) + Optikdetkl(:,:,end,:)*(Along(i)>A(end));
 HLdetkla(:,:,i,:)   = HLdetkla(:,:,i,:)*(Along(i)>A(1))*(Along(i)<A(end)) + Opthldetkl(:,:,1,:)*(Along(i)<A(1)) + Opthldetkl(:,:,end,:)*(Along(i)>A(end));
end


Kdetkla  = (1 + IKdetkla - repmat(permute(dK,[2 3 4 1]),[size(Kdet,2) size(Ldet,2) nA+extraA 1])).*repmat(permute(Kdeta,[2 4 1 3]),[1 size(Ldet,2) 1 1]);
Ldetkla  = (1 + HLdetkla - repmat(permute(dL,[2 3 4 1]),[size(Kdet,2) size(Ldet,2) nA+extraA 1])).*repmat(permute(Ldeta,[4 2 1 3]),[size(Kdet,2) 1 1 1]);

argdetK = nan(size(Kdet,2),size(Ldet,2),nA+extraA,nF);
argdetL = nan(size(Kdet,2),size(Ldet,2),nA+extraA,nF);
for a = 1 : nA + extraA
 [~,argdetK(:,:,a,:)] = min(abs(repmat(Kdetkla(:,:,a,:),[1 1 1 1 size(Kdet,2)]) - repmat(permute(Kdeta(a,:,:),[4 5 1 3 2]),[size(Kdet,2) size(Ldet,2) 1 1 1])),[],5);
 [~,argdetL(:,:,a,:)] = min(abs(repmat(Ldetkla(:,:,a,:),[1 1 1 1 size(Ldet,2)]) - repmat(permute(Ldeta(a,:,:),[4 5 1 3 2]),[size(Kdet,2) size(Ldet,2) 1 1 1])),[],5);
end 
Vdet = Vdetkla;

% Prepare simulation by creating some zero matrices.
 K_sim   = zeros(N,T,Kappa); L_sim   = zeros(N,T,Kappa);   K_elem  = zeros(N,T,Kappa);   L_elem  = zeros(N,T,Kappa); IK_sim = zeros(N,T,Kappa); 
 HL_sim  = zeros(N,T,Kappa); type    = zeros(N,T+1,Kappa); typesim = zeros(N,T+1,Kappa); Asim    = zeros(N,T+1,Kappa); Aord = zeros(N,T+1,Kappa); 
 V_sim   = zeros(N,T,Kappa); A_sim   = zeros(N,T+1,Kappa); dK_sim  = zeros(N,T,Kappa);   dL_sim  = zeros(N,T,Kappa); H_use = zeros(N,Tlast,Kappa); K_use1 = zeros(N,Tlast+1);
 IK_use  = zeros(N,Tlast,Kappa); HL_use = zeros(N,Tlast,Kappa); Y_use = zeros(N,Tlast,Kappa); V_use = zeros(N,Tlast,Kappa); I_use = zeros(N,Tlast,Kappa); L_use1 = zeros(N,Tlast+1);
 A_use   = zeros(N,Tlast,Kappa); K_use  = zeros(N,Tlast,Kappa); L_use = zeros(N,Tlast,Kappa); Y_gro = zeros(N,Tlast-1,Kappa); epsK = zeros(N,Tlast+1,Kappa); epsL = zeros(N,Tlast+1,Kappa);
 meanIK = zeros(Kappa,1); meanHL = zeros(Kappa,1); stdIK = zeros(Kappa,1); stdHL = zeros(Kappa,1); corrIK_HL = zeros(Kappa,1); 
 stdYgr = zeros(Kappa,1); corrIK_lagIK = zeros(Kappa,1); corrIK_lagHL = zeros(Kappa,1); corrHL_lagIK = zeros(Kappa,1);
 corrHL_lagHL   = zeros(Kappa,1); corrYgr_lagIK  = zeros(Kappa,1); corrYgr_lagHL   = zeros(Kappa,1); corrIK_lagYgr   = zeros(Kappa,1);
 corrHL_lagYgr  = zeros(Kappa,1); corrYgr_lagYgr = zeros(Kappa,1); corrIK_lag2IK   = zeros(Kappa,1); corrIK_lag2HL   = zeros(Kappa,1); 
 corrHL_lag2IK  = zeros(Kappa,1); corrHL_lag2HL  = zeros(Kappa,1); corrYgr_lag2IK  = zeros(Kappa,1); corrYgr_lag2HL  = zeros(Kappa,1); 
 corrIK_lag2Ygr = zeros(Kappa,1); corrHL_lag2Ygr = zeros(Kappa,1); corrYgr_lag2Ygr = zeros(Kappa,1); corrIK_Ygr      = zeros(Kappa,1); corrHL_Ygr = zeros(Kappa,1);
 IKneg = zeros(Kappa,1); IKzero = zeros(Kappa,1); IKpos = zeros(Kappa,1); HLneg = zeros(Kappa,1); HLzero = zeros(Kappa,1); HLpos = zeros(Kappa,1);
 dLneg = zeros(Kappa,1); dLzero = zeros(Kappa,1); dLpos = zeros(Kappa,1); Rnng  = zeros(Kappa,1); Rnzg   = zeros(Kappa,1); Rnpg  = zeros(Kappa,1); 
 Rzng  = zeros(Kappa,1); Rzzg   = zeros(Kappa,1); Rzpg  = zeros(Kappa,1); Rpng  = zeros(Kappa,1); Rpzg   = zeros(Kappa,1); Rppg  = zeros(Kappa,1);
 Rnnn  = zeros(Kappa,1); Rnzn   = zeros(Kappa,1); Rnpn  = zeros(Kappa,1); Rznn  = zeros(Kappa,1); Rzzn   = zeros(Kappa,1); Rzpn  = zeros(Kappa,1);   
 Rpnn  = zeros(Kappa,1); Rpzn   = zeros(Kappa,1); Rppn  = zeros(Kappa,1); ARnng = zeros(Kappa,1); ARnzg  = zeros(Kappa,1); ARzng = zeros(Kappa,1);  
 ARzzg = zeros(Kappa,1); ARnnn  = zeros(Kappa,1); ARnzn = zeros(Kappa,1); ARznn = zeros(Kappa,1); ARzzn  = zeros(Kappa,1);
 for s = 1 : Kappa;
  % Use different seed for each simulation round.
  seed = seed + s;
  % Draw uniformally distributed firm types and divide into nF types:
  if nF > 1;
    RandStream.setDefaultStream(RandStream('mt19937ar','seed',1*(27042011+s)));
    typesim(:,1,s) = rand(N,1);
    if nF == 2;
      type(:,1,s) = 1*(typesim(:,1,s)<0.5) + 2*(typesim(:,1,s)>=0.5);
    else
      type(:,1,s) = 1*(typesim(:,1,s)<1/3) + 2*(typesim(:,1,s)>=1/3)*(typesim(:,1,s)<2/3) + 3*(typesim(:,1,s)>2/3);
    end
    typesim(:,1,s) = Epsdepr(type(:,1,s));
    % Draw normally distributed shocks for the firm type transition process.
    RandStream.setDefaultStream(RandStream('mt19937ar','seed',2*(27042011+s)));
    v = sigdK.*randn(N,T+1);
    % Build up panel of firm types in terms of heterogenuous depr. rates.
    for t=2:T+1;
      typesim(:,t,s) = mean(dK).*(1-ro) + ro.*typesim(:,t-1,s) + v(:,t);
    end
    % These typesim()-components may be put into nF different categories.
    type(:,2:T+1,s) = (typesim(:,2:T+1,s)<Epsdepr(2));
    for j=2:nF-1;
      type(:,2:T+1,s) = type(:,2:T+1,s) + j.*(typesim(:,2:T+1,s)>=Epsdepr(j)).*(typesim(:,2:T+1,s)<=Epsdepr(j+1));
    end
    type(:,2:T+1,s) = type(:,2:T+1,s) + nF.*(typesim(:,2:T+1,s)>=Epsdepr(nF));
    else
    type(:,:,s) = ones(N,T+1,1);
  end
  
  % Draw normally distributed shocks for the firm level AR(1)-process.
  RandStream.setDefaultStream(RandStream('mt19937ar','seed',3*(27042011+s)));
  u = sigma.*randn(N,T+1);
  Asim(:,1,s) = u(:,1);
  % Build up panel of firm level productivity components.
  for t=2:T+1;
   Asim(:,t,s) = mu.*(1-rho) + rho.* Asim(:,t-1,s) + u(:,t);
  end
  % These Afirm_sim()-components may be put into nA+extraA different categories. 
  Aord(:,:,s) = (Asim(:,:,s)<Epsfirmlong(2));
  for j=2:nA+extraA-1;
    Aord(:,:,s) = Aord(:,:,s) + j.*(Asim(:,:,s)>=Epsfirmlong(j)).*(Asim(:,:,s)<=Epsfirmlong(j+1));
  end
  Aord(:,:,s) = Aord(:,:,s) + (nA+extraA).*(Asim(:,:,s)>Epsfirmlong(nA+extraA));
  
  
  % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
  % % %   Determine year by year investment rates and input levels    % % %
  % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
   
  IKvect   = reshape(IKdetkla,size(Kdet,2)*size(Ldet,2)*(nA+extraA)*nF,1,1,1);
  HLvect   = reshape(HLdetkla,size(Kdet,2)*size(Ldet,2)*(nA+extraA)*nF,1,1,1);
  Optkvect = reshape(argdetK,size(Kdet,2)*size(Ldet,2)*(nA+extraA)*nF,1,1,1);
  Optlvect = reshape(argdetL,size(Kdet,2)*size(Ldet,2)*(nA+extraA)*nF,1,1,1);
  Kdetvect = reshape(repmat(permute(Kdeta,[2 4 3 1]),[1 size(Ldet,2) 1 1]),size(Kdet,2)*size(Ldet,2)*(nA+extraA)*nF,1,1,1);
  Ldetvect = reshape(repmat(permute(Ldeta,[4 2 3 1]),[size(Kdet,2) 1 1 1]),size(Kdet,2)*size(Ldet,2)*(nA+extraA)*nF,1,1,1);
  Vvect    = reshape(Vdet,length(Kdet)*length(Ldet)*(nA+extraA)*nF,1,1,1);
  K_elem(:,1,s) = floor(length(Kdet)/2);
  L_elem(:,1,s) = floor(length(Ldet)/2);
  K_sim(:,1,s) = diag(Kdeta(Aord(:,1,s),K_elem(:,1,s)));
  L_sim(:,1,s) = diag(Ldeta(Aord(:,1,s),L_elem(:,1,s)));
  for t = 1 : T;
   A_sim(:,t,s)    = Along(Aord(:,t,s)');
   dK_sim(:,t,s)   = dK(type(:,t,s)');
   dL_sim(:,t,s)   = dL(type(:,t,s)');
   K_elem(:,t+1,s) = Optkvect((type(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2)*(nA+extraA) + (Aord(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2) + (L_elem(:,t,s)'-1)*size(Kdet,2) + K_elem(:,t,s)');
   L_elem(:,t+1,s) = Optlvect((type(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2)*(nA+extraA) + (Aord(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2) + (L_elem(:,t,s)'-1)*size(Kdet,2) + K_elem(:,t,s)');
   V_sim(:,t+1,s)  = Vvect((type(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2)*(nA+extraA) + (Aord(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2) + (L_elem(:,t,s)'-1)*size(Kdet,2) + K_elem(:,t,s)');
   K_sim(:,t+1,s)  = Kdetvect((type(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2)*(nA+extraA) + (Aord(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2) + (L_elem(:,t+1,s)'-1)*size(Kdet,2) + K_elem(:,t+1,s)');
   L_sim(:,t+1,s)  = Ldetvect((type(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2)*(nA+extraA) + (Aord(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2) + (L_elem(:,t+1,s)'-1)*size(Kdet,2) + K_elem(:,t+1,s)');
   IK_sim(:,t,s)   = IKvect((type(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2)*(nA+extraA) + (Aord(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2) + (L_elem(:,t,s)'-1)*size(Kdet,2) + K_elem(:,t,s)');
   HL_sim(:,t,s)   = HLvect((type(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2)*(nA+extraA) + (Aord(:,t,s)'-1)*size(Kdet,2)*size(Ldet,2) + (L_elem(:,t,s)'-1)*size(Kdet,2) + K_elem(:,t,s)');   
   K_elem(:,t+1,s) = max(K_elem(:,t+1,s) - round(2*(detail+1)*aK*((log(Along(Aord(:,t+1,s))) - log(Along(Aord(:,t,s))))/(log(A(end))-log(A(1))))) - (detail+1)*fK*(type(:,t+1,s) - type(:,t,s)),1);
   L_elem(:,t+1,s) = max(L_elem(:,t+1,s) - round(2*(detail+1)*aL*((log(Along(Aord(:,t+1,s))) - log(Along(Aord(:,t,s))))/(log(A(end))-log(A(1))))) - (detail+1)*fK*(type(:,t+1,s) - type(:,t,s)),1);
   % K_elem(:,t+1,s) = K_elem(:,t+1,s) - round(2*(detail+1)*aK*((log(Along(Aord(:,t+1,s))) - log(Along(Aord(:,t,s))))/(log(A(end))-log(A(1))))) - (detail+1)*fK*(type(:,t+1,s) - type(:,t,s));
   % L_elem(:,t+1,s) = L_elem(:,t+1,s) - round(2*(detail+1)*aL*((log(Along(Aord(:,t+1,s))) - log(Along(Aord(:,t,s))))/(log(A(end))-log(A(1))))) - (detail+1)*fK*(type(:,t+1,s) - type(:,t,s));
   K_sim(:,t+1,s)  = diag(Kdeta(Aord(:,t+1,s),K_elem(:,t+1,s)));
   L_sim(:,t+1,s)  = diag(Ldeta(Aord(:,t+1,s),L_elem(:,t+1,s)));
  end

 
  % Keep the T_last periods of the simulated data and add measurement errors:   
  RandStream.setDefaultStream(RandStream('mt19937ar','seed',4*(27042011+s)));
  epsIK(:,:,s)   = me(1).*randn(N,Tlast);
  RandStream.setDefaultStream(RandStream('mt19937ar','seed',5*(27042011+s)));
  epsHL(:,:,s)   = me(2).*randn(N,Tlast);
  A_use(:,:,s)  = A_sim(:,(T - Tlast + 1): T,s);
  K_use(:,:,s)  = K_sim(:,(T-Tlast+1):T,s);
  L_use(:,:,s)  = L_sim(:,(T-Tlast+1):T,s);
  Y_use(:,:,s)  = A_use(:,:,s).*K_use(:,:,s).^etaK.*L_use(:,:,s).^etaL;
  IK_use(:,:,s) = IK_sim(:,(T-Tlast+1):T,s).*exp(epsIK(:,:,s));
  HL_use(:,:,s) = HL_sim(:,(T-Tlast+1):T,s).*exp(epsHL(:,:,s));
  V_use(:,:,s)  = V_sim(:,(T - Tlast + 1): T,s);
  Y_gro(:,:,s)  = (Y_use(:,2:end,s)./Y_use(:,1:end-1,s) - 1);
  YK_use(:,:,s) = Y_use(:,:,s)./K_use(:,:,s);
  YL_use(:,:,s) = Y_use(:,:,s)./L_use(:,:,s);
  dK_use(:,:,s) = dK_sim(:,(T-Tlast+1):T,s);
  dL_use(:,:,s) = dL_sim(:,(T-Tlast+1):T,s);
  
  
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% % % % %       Calculate moments in simulated samples.       % % % % % %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
    
  % Numbered list of moments is given in getmoments.m.

  % Determine some inaction definitions:
  IKinaction = [-0.03 0.03];
  HLgrossinaction = [-0.03 0.03];
  HLnetinaction = [0.07 0.13];
  
  % Record some moments:
  meanIK(s,1) = mean(reshape(IK_use(:,:,s),size(IK_use,1)*size(IK_use,2),1));
  stdIK(s,1)  = std(reshape(IK_use(:,:,s),size(IK_use,1)*size(IK_use,2),1));
  meanHL(s,1) = mean(reshape(HL_use(:,:,s),size(IK_use,1)*size(IK_use,2),1));
  stdHL(s,1)  = std(reshape(HL_use(:,:,s),size(IK_use,1)*size(IK_use,2),1));
  stdYgr(s,1) = std(reshape(Y_gro(:,:,s),size(Y_gro,1)*size(Y_gro,2),1));
  skewIK(s,1) = skewness(reshape(IK_use(:,:,s),size(IK_use,1)*size(IK_use,2),1));
  skewHL(s,1) = skewness(reshape(HL_use(:,:,s),size(HL_use,1)*size(HL_use,2),1));
  IKneg(s,1)  = mean(reshape(IK_use(:,:,s)<IKinaction(1),size(IK_use,1)*size(IK_use,2),1));
  IKzero(s,1) = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  IKpos(s,1)  = mean(reshape(IK_use(:,:,s)>IKinaction(2),size(IK_use,1)*size(IK_use,2),1));
  HLneg(s,1)  = mean(reshape(HL_use(:,:,s)<HLgrossinaction(1),size(IK_use,1)*size(IK_use,2),1));
  HLzero(s,1) = mean(reshape((HL_use(:,:,s)>=HLgrossinaction(1)).*(HL_use(:,:,s)<=HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  HLpos(s,1)  = mean(reshape(HL_use(:,:,s)>HLgrossinaction(2),size(IK_use,1)*size(IK_use,2),1));
  dLneg(s,1)  = mean(reshape(HL_use(:,:,s)<HLnetinaction(1),size(IK_use,1)*size(IK_use,2),1));
  dLzero(s,1) = mean(reshape((HL_use(:,:,s)>=HLnetinaction(1)).*(HL_use(:,:,s)<=HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  dLpos(s,1)  = mean(reshape(HL_use(:,:,s)>HLnetinaction(2),size(IK_use,1)*size(IK_use,2),1));
  Rnng(s,1)   = mean(reshape((IK_use(:,:,s)<IKinaction(1)).*(HL_use(:,:,s)<HLgrossinaction(1)),size(IK_use,1)*size(IK_use,2),1));
  Rnzg(s,1)   = mean(reshape((IK_use(:,:,s)<IKinaction(1)).*(HL_use(:,:,s)>=HLgrossinaction(1)).*(HL_use(:,:,s)<=HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rnpg(s,1)   = mean(reshape((IK_use(:,:,s)<IKinaction(1)).*(HL_use(:,:,s)>HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rzng(s,1)   = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)<HLgrossinaction(1)),size(IK_use,1)*size(IK_use,2),1));
  Rzzg(s,1)   = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)>=HLgrossinaction(1)).*(HL_use(:,:,s)<=HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rzpg(s,1)   = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)>HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rpng(s,1)   = mean(reshape((IK_use(:,:,s)>IKinaction(2)).*(HL_use(:,:,s)<HLgrossinaction(1)),size(IK_use,1)*size(IK_use,2),1));
  Rpzg(s,1)   = mean(reshape((IK_use(:,:,s)>IKinaction(2)).*(HL_use(:,:,s)>=HLgrossinaction(1)).*(HL_use(:,:,s)<=HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rppg(s,1)   = mean(reshape((IK_use(:,:,s)>IKinaction(2)).*(HL_use(:,:,s)>HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rnnn(s,1)   = mean(reshape((IK_use(:,:,s)<IKinaction(1)).*(HL_use(:,:,s)<HLnetinaction(1)),size(IK_use,1)*size(IK_use,2),1)); 
  Rnzn(s,1)   = mean(reshape((IK_use(:,:,s)<IKinaction(1)).*(HL_use(:,:,s)>=HLnetinaction(1)).*(HL_use(:,:,s)<=HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rnpn(s,1)   = mean(reshape((IK_use(:,:,s)<IKinaction(1)).*(HL_use(:,:,s)>HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rznn(s,1)   = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)<HLnetinaction(1)),size(IK_use,1)*size(IK_use,2),1));
  Rzzn(s,1)   = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)>=HLnetinaction(1)).*(HL_use(:,:,s)<=HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rzpn(s,1)   = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)>HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rpnn(s,1)   = mean(reshape((IK_use(:,:,s)>IKinaction(2)).*(HL_use(:,:,s)<HLnetinaction(1)),size(IK_use,1)*size(IK_use,2),1));
  Rpzn(s,1)   = mean(reshape((IK_use(:,:,s)>IKinaction(2)).*(HL_use(:,:,s)>=HLnetinaction(1)).*(HL_use(:,:,s)<=HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  Rppn(s,1)   = mean(reshape((IK_use(:,:,s)>IKinaction(2)).*(HL_use(:,:,s)>HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  ARnng(s,1)  = mean(reshape(((IK_use(:,:,s)<IKinaction(1))+(IK_use(:,:,s)>IKinaction(2))).*((HL_use(:,:,s)<HLgrossinaction(1))+(HL_use(:,:,s)>HLgrossinaction(2))),size(IK_use,1)*size(IK_use,2),1));
  ARnzg(s,1)  = mean(reshape(((IK_use(:,:,s)<IKinaction(1))+(IK_use(:,:,s)>IKinaction(2))).*(HL_use(:,:,s)>=HLgrossinaction(1)).*(HL_use(:,:,s)<=HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  ARzng(s,1)  = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*((HL_use(:,:,s)<HLgrossinaction(1))+(HL_use(:,:,s)>HLgrossinaction(2))),size(IK_use,1)*size(IK_use,2),1));
  ARzzg(s,1)  = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)>=HLgrossinaction(1)).*(HL_use(:,:,s)<=HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  ARnnn(s,1)  = mean(reshape(((IK_use(:,:,s)<IKinaction(1))+(IK_use(:,:,s)>IKinaction(2))).*((HL_use(:,:,s)<HLnetinaction(1))+(HL_use(:,:,s)>HLnetinaction(2))),size(IK_use,1)*size(IK_use,2),1));
  ARnzn(s,1)  = mean(reshape(((IK_use(:,:,s)<IKinaction(1))+(IK_use(:,:,s)>IKinaction(2))).*(HL_use(:,:,s)>=HLnetinaction(1)).*(HL_use(:,:,s)<=HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  ARznn(s,1)  = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*((HL_use(:,:,s)<HLnetinaction(1))+(HL_use(:,:,s)>HLnetinaction(2))),size(IK_use,1)*size(IK_use,2),1));
  ARzzn(s,1)  = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)>=HLnetinaction(1)).*(HL_use(:,:,s)<=HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  
  % Determine some inaction definitions:
  IKinaction = [-0.02 0.02];
  HLgrossinaction = [-0.02 0.02];
  HLnetinaction = [0.08 0.12];

  IKneg2(s,1)       = mean(reshape(IK_use(:,:,s)<IKinaction(1),size(IK_use,1)*size(IK_use,2),1));
  IKzero2(s,1)      = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  HLzero2(s,1)      = mean(reshape((HL_use(:,:,s)>=HLgrossinaction(1)).*(HL_use(:,:,s)<=HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  dLzero2(s,1)      = mean(reshape((HL_use(:,:,s)>=HLnetinaction(1)).*(HL_use(:,:,s)<=HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  simzerogross(s,1) = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)>=HLgrossinaction(1)).*(HL_use(:,:,s)<=HLgrossinaction(2)),size(IK_use,1)*size(IK_use,2),1));
  simzeronet(s,1)   = mean(reshape((IK_use(:,:,s)>=IKinaction(1)).*(IK_use(:,:,s)<=IKinaction(2)).*(HL_use(:,:,s)>=HLnetinaction(1)).*(HL_use(:,:,s)<=HLnetinaction(2)),size(IK_use,1)*size(IK_use,2),1));

  C = crosscorr(reshape(IK_use(:,:,s)',N*Tlast,1),reshape(HL_use(:,:,s)',N*Tlast,1),2,[]);
  corrIK_HL(s,1)     = C(3);
  corrIK_lagHL(s,1)  = C(2);
  corrIK_lag2HL(s,1) = C(1);
  corrHL_lagIK(s,1)  = C(4);
  corrHL_lag2IK(s,1) = C(5);
  C = crosscorr(reshape(IK_use(:,2:end,s)',N*(Tlast-1),1),reshape(Y_gro(:,:,s)',N*(Tlast-1),1),2,[]);
  corrIK_Ygr(s,1)     = C(3);
  corrIK_lagYgr(s,1)  = C(2);
  corrIK_lag2Ygr(s,1) = C(1);
  corrYgr_lagIK(s,1)  = C(4);
  corrYgr_lag2IK(s,1) = C(5);
  C = crosscorr(reshape(HL_use(:,2:end,s)',N*(Tlast-1),1),reshape(Y_gro(:,:,s)',N*(Tlast-1),1),2,[]);
  corrHL_Ygr(s,1)     = C(3);
  corrHL_lagYgr(s,1)  = C(2);
  corrHL_lag2Ygr(s,1) = C(1);
  corrYgr_lagHL(s,1)  = C(4);
  corrYgr_lag2HL(s,1) = C(5);
  C = crosscorr(reshape(K_use(:,1:end,s)',N*Tlast,1),reshape(Y_use(:,:,s)',N*Tlast,1),2,[]);
  corrK_Y(s,1)        = C(3);
  C = crosscorr(reshape(L_use(:,1:end,s)',N*Tlast,1),reshape(Y_use(:,:,s)',N*Tlast,1),2,[]);
  corrL_Y(s,1)        = C(3);
  C = crosscorr(reshape(K_use(:,1:end,s)',N*Tlast,1),reshape(L_use(:,:,s)',N*Tlast,1),2,[]);
  corrK_L(s,1)        = C(3);
  C = crosscorr(reshape(IK_use(:,1:end,s)',N*Tlast,1),reshape(log(YK_use(:,:,s))',N*Tlast,1),2,[]);
  corrIK_YK(s,1)      = C(3);
  corrIK_lagYK(s,1)   = C(2);
  corrIK_lag2YK(s,1)  = C(1);
  C = crosscorr(reshape(HL_use(:,1:end,s)',N*Tlast,1),reshape(log(YL_use(:,:,s))',N*Tlast,1),2,[]);
  corrHL_YL(s,1)      = C(3);
  corrHL_lagYL(s,1)   = C(2);
  corrHL_lag2YL(s,1)  = C(1);
  C = crosscorr(reshape(HL_use(:,1:end,s)',N*Tlast,1),reshape(log(YK_use(:,:,s))',N*Tlast,1),2,[]);
  corrHL_YK(s,1)      = C(3);
  corrHL_lagYK(s,1)   = C(2);
  corrHL_lag2YK(s,1)  = C(1);
  C = crosscorr(reshape(Y_gro(:,1:end,s)',N*(Tlast-1),1),reshape(log(YK_use(:,2:end,s))',N*(Tlast-1),1),2,[]);
  corrYgr_YK(s,1)     = C(3);
  corrYgr_lagYK(s,1)  = C(2);
  corrYgr_lag2YK(s,1) = C(1);
  C = crosscorr(reshape(Y_gro(:,1:end,s)',N*(Tlast-1),1),reshape(log(YL_use(:,2:end,s))',N*(Tlast-1),1),2,[]);
  corrYgr_YL(s,1)     = C(3);
  corrYgr_lagYL(s,1)  = C(2);
  corrYgr_lag2YL(s,1) = C(1);
  C = autocorr(reshape(IK_use(:,:,s)',N*Tlast,1),2,[]);
  corrIK_lagIK(s,1)  = C(2);
  corrIK_lag2IK(s,1) = C(3);
  C = autocorr(reshape(HL_use(:,:,s)',N*Tlast,1),2,[]);
  corrHL_lagHL(s,1)  = C(2);
  corrHL_lag2HL(s,1) = C(3);
  C = autocorr(reshape(Y_gro(:,:,s)',N*(Tlast-1),1),2,[]);
  corrYgr_lagYgr(s,1)  = C(2);
  corrYgr_lag2Ygr(s,1) = C(3);
  C = autocorr(reshape(K_use(:,:,s)',N*Tlast,1),2,[]);
  corrK_lagK(s,1) = C(2);
  corrK_lag2K(s,1) = C(3);
  C = autocorr(reshape(L_use(:,:,s)',N*Tlast,1),2,[]);
  corrL_lagL(s,1) = C(2);
  corrL_lag2L(s,1) = C(3);
  C = autocorr(reshape(log(YK_use(:,:,s))',N*Tlast,1),2,[]);
  corrYK_lagYK(s,1) = C(2);
  corrYK_lag2YK(s,1) = C(3);
  C = autocorr(reshape(log(YL_use(:,:,s))',N*Tlast,1),2,[]);
  corrYL_lagYL(s,1) = C(2);
  corrYL_lag2YL(s,1) = C(3);
  
  
  end  % END LOOP OVER SIMULATION ROUNDS 1 : KAPPA.   
  
 % Create moment matrix with dimensions {(Kappa) x (# of moments)}.
  mom = [corrIK_lagIK  corrIK_lagHL  corrIK_lagYgr  corrHL_lagIK  corrHL_lagHL  corrHL_lagYgr  corrYgr_lagIK  corrYgr_lagHL  corrYgr_lagYgr  ...
         corrIK_lag2IK corrIK_lag2HL corrIK_lag2Ygr corrHL_lag2IK corrHL_lag2HL corrHL_lag2Ygr corrYgr_lag2IK corrYgr_lag2HL corrYgr_lag2Ygr ...
         corrIK_HL     corrIK_Ygr    corrHL_Ygr     meanIK        stdIK         meanHL         stdHL          stdYgr         IKneg           ...
         IKzero        IKpos         HLneg          HLzero        HLpos         dLneg          dLzero         dLpos          Rnng            ...
         Rnzg          Rnpg          Rzng           Rzzg          Rzpg          Rpng           Rpzg           Rppg           Rnnn            ...
         Rnzn          Rnpn          Rznn           Rzzn          Rzpn          Rpnn           Rpzn           Rppn           ARnng           ...
         ARnzg         ARzng         ARzzg          ARnnn         ARnzn         ARznn          ARzzn          corrK_lagK     corrK_lag2K     ...
         corrL_lagL    corrL_lag2L   corrK_Y        corrL_Y       corrK_L       skewIK         corrIK_YK      corrIK_lagYK   corrIK_lag2YK   ...
         corrHL_YL     corrHL_lagYL  corrHL_lag2YL  corrHL_YK     corrHL_lagYK  corrHL_lag2YK  corrYgr_YK     corrYgr_lagYK  corrYgr_lag2YK  ...
         corrYgr_YL    corrYgr_lagYL corrYgr_lag2YL corrYK_lagYK  corrYK_lag2YK corrYL_lagYL   corrYL_lag2YL  skewHL         IKneg2          ...
         IKzero2       HLzero2       dLzero2        simzerogross  simzeronet];     
    
    
 % Return moment vector as the means over Kappa rounds of simulations.
 simmom = mean(mom);
 
 % Compare with empirical moments already obtained and loaded into memory.
 criterion = evaluation(empmom,simmom,covar,momlist); 
 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%             Display results if execution mode = 'Sim'               %
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

if strcmp(mode,'Sim');
    fprintf('--------------------------------------\n')
    fprintf('         Simulation inputs \n')
    fprintf('--------------------------------------\n')
    fprintf('Cost parameters \n')
    fprintf(' bK         = %3.5f \n', params(1))
    fprintf(' bL         = %3.5f \n', params(2))
    fprintf(' alphaK     = %3.5f \n', params(3))
    fprintf(' alphaL     = %3.5f \n', params(4))
    fprintf(' alphaKL    = %3.5f \n', params(5))
    fprintf(' resaleloss = %3.5f \n', params(6))
    fprintf(' rho        = %3.5f \n', params(7))
    fprintf(' sigma      = %3.5f \n', params(8))
    fprintf('--------------------------------------\n')
    fprintf('Other parameters \n')
    fprintf(' n(K)       = %3.0f \n',nK)
    fprintf(' n(L)       = %3.0f \n',nL)
    fprintf(' n(A)       = %3.0f (%1.0f) \n',nA,nA+extraA)
    fprintf(' n(F)       = %3.0f \n',nF)
    fprintf(' min(dK)    = %3.4f \n',min(dK))
    fprintf(' max(dK)    = %3.4f \n',max(dK))
    fprintf(' meK        = %3.4f \n',params(12))
    fprintf(' meL        = %3.4f \n',params(13))
    fprintf('----------------------------------------------------\n')
    fprintf('           Modes \n')
    fprintf('----------------------------------------------------\n')
    if strcmp(Lmode,'LaborGross')
     fprintf('Fixed costs of labor associated with gross hirings \n')
    else
     fprintf('Fixed costs of labor associated with net hirings \n')
    end
    if strcmp(FC,'Fixed')
     fprintf('All fixed costs are purely fixed \n')
    else
     fprintf('All fixed costs are scaled by plant size \n')
    end
    fprintf('----------------------------------------------------\n')
    fprintf(' \n')
    fprintf('-----------------------------------------------------------\n')
    fprintf('                  Simulation results \n')
    fprintf('-----------------------------------------------------------\n')   
    fprintf('                     Empirical moments  Simulation moments \n')
    fprintf('Corr(IK(t),IK(t-1))       %7.4f           %7.4f \n',empmom(1),simmom(1))
    fprintf('Corr(IK(t),HL(t-1))       %7.4f           %7.4f \n',empmom(2),simmom(2))
    fprintf('Corr(IK(t),Ygr(t-1))      %7.4f           %7.4f \n',empmom(3),simmom(3))
    fprintf('Corr(HL(t),IK(t-1))       %7.4f           %7.4f \n',empmom(4),simmom(4))
    fprintf('Corr(HL(t),HL(t-1))       %7.4f           %7.4f \n',empmom(5),simmom(5))
    fprintf('Corr(HL(t),Ygr(t-1))      %7.4f           %7.4f \n',empmom(6),simmom(6))
    fprintf('Corr(Ygr(t),IK(t-1))      %7.4f           %7.4f \n',empmom(7),simmom(7))
    fprintf('Corr(Ygr(t),HL(t-1))      %7.4f           %7.4f \n',empmom(8),simmom(8))
    fprintf('Corr(Ygr(t),Ygr(t-1))     %7.4f           %7.4f \n',empmom(9),simmom(9))
    fprintf('Corr(IK(t),IK(t-2))       %7.4f           %7.4f \n',empmom(10),simmom(10))
    fprintf('Corr(IK(t),HL(t-2))       %7.4f           %7.4f \n',empmom(11),simmom(11))
    fprintf('Corr(IK(t),Ygr(t-2))      %7.4f           %7.4f \n',empmom(12),simmom(12))
    fprintf('Corr(HL(t),IK(t-2))       %7.4f           %7.4f \n',empmom(13),simmom(13))
    fprintf('Corr(HL(t),HL(t-2))       %7.4f           %7.4f \n',empmom(14),simmom(14))
    fprintf('Corr(HL(t),Ygr(t-2))      %7.4f           %7.4f \n',empmom(15),simmom(15))
    fprintf('Corr(Ygr(t),IK(t-2))      %7.4f           %7.4f \n',empmom(16),simmom(16))
    fprintf('Corr(Ygr(t),HL(t-2))      %7.4f           %7.4f \n',empmom(17),simmom(17))
    fprintf('Corr(Ygr(t),Ygr(t-2))     %7.4f           %7.4f \n',empmom(18),simmom(18))
    fprintf('Std(IK)                   %7.4f           %7.4f \n',empmom(23),simmom(23))
    fprintf('Std(HL)                   %7.4f           %7.4f \n',empmom(25),simmom(25))
    fprintf('Std(Ygr)                  %7.4f           %7.4f \n',empmom(26),simmom(26))
    fprintf('Corr(K(t),K(t-1)          %7.4f           %7.4f \n',empmom(62),simmom(62))
    fprintf('Corr(K(t),K(t-2)          %7.4f           %7.4f \n',empmom(63),simmom(63))
    fprintf('Corr(L(t),L(t-1)          %7.4f           %7.4f \n',empmom(64),simmom(64))
    fprintf('Corr(L(t),L(t-2)          %7.4f           %7.4f \n',empmom(65),simmom(65))
    fprintf('Corr(lnK(t),lnY(t)        %7.4f           %7.4f \n',empmom(66),simmom(66))
    fprintf('Corr(lnL(t),lnY(t)        %7.4f           %7.4f \n',empmom(67),simmom(67))
    fprintf('Corr(lnK(t),lnL(t)        %7.4f           %7.4f \n',empmom(68),simmom(68))
    fprintf('-----------------------------------------------------------\n')   
    fprintf('\n')
    fprintf('----------------------------- \n')
    fprintf('Criterion = %7.1f \n', criterion)
    fprintf('------------------------------\n')   
    fprintf('\n')
end

end

end
